MODULE md1_LIP_coeffs

   USE md1_global_consts
   USE md1_global_types
   IMPLICIT NONE
   PRIVATE
   PUBLIC :: md1_init_LIP_coeffs,            &
             md1_free_LIP_coeffs,            &
             md1_get_LIP_coeffs

   REAL(REALK), PARAMETER :: LIP1_coeffs(4) = (/   &
                ! f0
                             one,                  &
                             -one,                 &
                ! f1
                             zero,                 &
                             one              /)
                             
   REAL(REALK), PARAMETER :: LIP2_coeffs(9) = (/   &
                ! f0
                             one,                  &
                             -three/two,           &
                             half,                 &
                ! f1
                             zero,                 &
                             two,                  &
                             -one,                 &
                ! f2
                             zero,                 &
                             -half,                &
                             half             /)

   REAL(REALK), PARAMETER :: LIP3_coeffs(16) = (/   &
                ! f0
                             one,                  &
                             -11d0/6d0,            &
                             one,                  &
                             -one/6d0,             &
                ! f1
                             zero,                 &
                             three,                &
                             -5d0/2d0,             &
                             half,                 &
                ! f2
                             zero,                 &
                             -3d0/2d0,             &
                             two,                  &
                             -half,                &
                ! f3
                             zero,                 &
                             one/three,            &
                             -half,                &
                             one/6d0          /)

   REAL(REALK), ALLOCATABLE, SAVE, TARGET :: LIP_coeffs(:,:,:)

   CONTAINS

!------------------------------------------------------------------------------
   
   SUBROUTINE md1_init_LIP_coeffs(order,h)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN) :: order
      REAL(REALK),   INTENT(IN) :: h(3)

      INTEGER(INTK) :: fn, xn, n, xyz
      REAL(REALK)   :: tmp(0:order,0:order)

      IF (ALLOCATED (LIP_coeffs)) DEALLOCATE (LIP_coeffs)
      ALLOCATE (LIP_coeffs(3,0:order,0:order))

      n = order + 1
      SELECT CASE (order)
        CASE (1)
          tmp(0:order,0:order) = RESHAPE(LIP1_coeffs,(/n,n/))
        CASE (2)
          tmp(0:order,0:order) = RESHAPE(LIP2_coeffs,(/n,n/))
        CASE (3)
          tmp(0:order,0:order) = RESHAPE(LIP3_coeffs,(/n,n/))
        CASE DEFAULT 
          STOP 'fixme'
      END SELECT

      ! Rescale according to h-value (grid spacing)
      DO xyz = 1,3
         DO fn = 0, order
            DO xn = 0, order
               LIP_coeffs(xyz,xn,fn) = tmp(xn,fn) / (h(xyz)**xn)
            END DO
         END DO
         LIP_coeffs(xyz,:,:) = TRANSPOSE (LIP_coeffs(xyz,:,:))
      END DO
   
   END SUBROUTINE md1_init_LIP_coeffs

!------------------------------------------------------------------------------

   SUBROUTINE md1_free_LIP_coeffs
      IMPLICIT NONE
      IF (ALLOCATED (LIP_coeffs)) DEALLOCATE (LIP_coeffs)
   END SUBROUTINE md1_free_LIP_coeffs

!------------------------------------------------------------------------------

   SUBROUTINE md1_get_LIP_coeffs(xyz,order,coeffs) 

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN) :: xyz, order
      REAL(REALK),   POINTER    :: coeffs(:,:)

      coeffs => LIP_coeffs(xyz,0:order,0:order)

   END SUBROUTINE md1_get_LIP_coeffs

!------------------------------------------------------------------------------

END MODULE md1_LIP_coeffs

